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We theoretically investigate light propagation and Anderson localization in one-dimensional disordered su¬ 
perlattices composed of dielectric stacks with graphene sheets in between. Disorder is introduced either on 
graphene material parameters {e.g. Fermi energy) or on the widths of the dielectric stacks. We derive an an¬ 
alytic expression for the localization length and compare it to numerical simulations using transfer matrix 
technique; a very good agreement is found. We demonstrate that the presence of graphene may strongly attenu¬ 
ate the anomalously delocalised Breswter modes, and is at the origin of a periodic dependence of ^ on frequency, 
in contrast to the usual asymptotic decay, oc By unveiling the effects of graphene on Anderson localiza¬ 
tion of light, we pave the way for new applications of graphene-based, disordered photonic devices in the THz 
spectral range. 

PACS numbers: 


I. INTRODUCTION 

Due to its extraordinary electronic and optical properties, 
graphene has emerged as an alternative material platform for 
applications in photonics and optoelectronic A partial but 
by no means exhaustive list of applications of graphene in 
photonics include high-speedphotodetectorP, optical mod- 
ulatorP, plasmonic deviceJ^® and ultrafast laserP. In ad¬ 
dition, graphene is a promising candidate to overcome one 
of the major existing hurdles to bring optics and electron¬ 
ics together, namely the efficient conversion between optical 
and electronic signals. Indeed, this can be facilitated by the 
fact that graphene enables strong, electric field-tunable optical 
transitions, and resonantly enhances light-mater interactions 
in sub-wavelength volumes. In practice this can be achieved, 
for instance, by i nteg rating a graphene layer into a photonic 
crystal nanocavity^. The presence of graphene also allows 
for an efficient electro-optical modul ation of photonic crys¬ 
tals nanocavities by electrostatic gatin^HE^. However, the in¬ 
tegration of graphene into photonic crystals is naturally prone 
to unavoidable disorder associated to the fabrication process. 
This constitutes per se a motivation to investigate the effects 
of disorder in photonic crystals containing graphene layers 
which, as far as we know, have not been considered in the 
literature so far. In addition to this technological and practical 
motivation, there is a very fundamental one as well, namely to 
understand the impact of graphene on Anderson localization 
of light. 

The concept of Anderson localization (AL) was originally 
conceived in the realm of condensed matter physics as a dis¬ 
order driven metal-insulator transitiorPI. Being an interfer¬ 
ence wave phenomenon, this concept has been extended to 
lighP^, acoustic waveJI^, and even Bose-Einstein condensed 
matter wave^. As a result, Anderson localization is today a 
truly interdisciplinary topic, and important contributions have 
emerged from different areas, ranging from condensed mat¬ 


ter, photonics, acoustics, atomic physics, and seismolog}C^. 
Dimensionality is crucial to AL, and in ID the vast majority 
of states is exponentially localised on a length scale given by 
the localization length regardless of the disorder strength. 
In optical systems exceptions do exist, and delocalised modes 
may occur in low-dimensional systems as a result of the pres¬ 
ence of correlation^^, necklace modeJ^^, or metamaterials 
with negative refractiorPH^. The question of whether these 
anomalies occur when graphene is integrated into disordered 
optical superlattices remains an open question. 

Bearing in mind both these technological and fundamental 
motivations, in the present paper we undertake an analytical 
and numerical investigation of Anderson localization of light 
in one-dimensional disordered superlattices composed of di¬ 
electric stacks with graphene layers in between, as depicted 
in Fig. [2 We consider two possible, realistic ways to model 
disorder: compositional and structural disorder. In the for¬ 
mer case disorder is introduced in graphene’s material param¬ 
eters, such as the Fermi energy, whereas in the latter the di¬ 
electric components of the superlattice have random widths. 
In both cases, we derive an analytic expression for the lo¬ 
calization length and compare it to numerical simulations 
using a transfer matrix technique; an overall very good agree¬ 
ment is found. In the case where the medium impedances 
match, we find that ^ exhibits an oscillatory behaviour as a 
function of frequency cc, in contrast to the usual asymptotic 
decay ^ ex . We demonstrate that graphene may strongly 
suppress the anomalously delocalised Brewster modes, as it 
induces additional refiexions at the superlattice interfaces. We 
also investigate the effects of inter and intraband transitions of 
the graphene conductivity on identifying the regimes where 
Anderson localization and absorption dominates light trans¬ 
mission. 

This paper is organised as follows. In Sec. II we present the 
analytical results, where we derive an expression for the local¬ 
ization length of disordered superlattices containing graphene 
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sheets. In Sec. Ill we present and discuss the numerical sim¬ 
ulations, based on transfer matrix technique, which are also 
compared to the analytical calculations. Finally, Sec. IV is 
devoted to the concluding remarks. We also present a number 
of appendices giving the details of the calculations and aim¬ 
ing at making the the text as self-contained as possible . To 
our best knowledge, there are only two published paper^^^^ 
dealing with similar problems to the one we consider in this 
paper, but in the context a metals, in which case only Drude’s 
conductivity plays a role. 


graphene 
X y_X 



FIG. 1. Color on-line. Schematic representation of the system. 


II. ANALYTICAL CALCULATION OF THE 
LOCALIZATION LENGTH 

Light promgation in a ID superlattice containing graphene 
layers (Fig. [l]) is modelled by the transfer matrix formalisnP^. 
The } transfer matrix connects the fields at the 

right of the n— th unit cell to those at left according to: 

^n+l _ 

where (V^l) refers to the right 

(left) propagating field in the n— th cell. For transverse elec¬ 
tric (TE) and transvere magnetic (TM) modes, ip refers to the 
electric and magnetic field, respectively. We consider the par¬ 
ticular case where detM’^ = 1, which occurs for systems 
with preserved time reversal symmetiyE^. In this case, one 
can show that may be written as 

i\/rn _ ( cosh^ye*"^^ sinh^ye^'^s \ 

where (pf are parameters that depend on the composition of 
the n— th cell, (from here on we omit the n dependence in (pi, 
except when strictly necessary to avoid any confusion.) For 
periodic systems with preserved time-reversed symmetry, cpi 
are real numbers and all the M’^’s are equal. We thus write 
= M^. One can write the photonic dispersion relatiorPl 
as cosy = (m^i + 17122)I where 

cos 7 = cosh (pi cos (p2 . ( 3 ) 


Disorder is introduced in the parameters (pi : 

c^i=4>1 + Hi ( 4 ) 


where 5 (pi describes random fiuctuations around the average 
value, and which may have different origins, as it will be de¬ 
tailed later in the paper. For a periodic system, a transforma¬ 
tion Mtransf = M;ircieM-eai (sce appendix]^ exists that maps 
the variables ip^ ^ into a new set of variables, denoted by 

Qn and Pn, such that = [Q^ = Mransf['0S 

These n ew variabl es describe a circle in phase spac^, with 
radius proportional to the electric field amplitude. 

Applying this transformation to Eq. Q, the transformed ma- 
trix M' = reads 


M' = 


(En 

\Gn 



( 5 ) 


where: 


En = cosh (pi cos (p2 — sinh (pi sin S(ps , (6) 

Fn = —(cosh (pi sin (p2 + sinh (pi cos S(ps) , ( 7 ) 

Gn = (cosh (pi sin (p2 — sinh (pi cos S(ps ), (8) 

Hn = cosh (pi cos (p2 + sinh (pi sin S(ps . ( 9 ) 

with V and r defined in appendix When (pi = (p^, we 
have d(p2, = 0 and Eqs. ([^-([^ lead to = cosy 

and En = — Gn = siny. When weak disorder is intro¬ 
duced, the trajectory of the points {Qn^Pn) results in a per¬ 
turbation of the circle. The recurrence equations defined by 

Xn 7 l ^ 

are similar to a Hamiltonian map of the clas¬ 
sical harmonic oscillator subjected to a parametric impulsive 
forc^, where Qn and Pn are the coordinate and conjugated 
moments, respectively, and y is the phase between successive 
kicks. 

The presence of disorder introduces a key length scale, the 
localization length In ID electronic systems all eigenmodes 
are exponentially localised, altho^h some exceptions do ex¬ 
ist in the realm of optical systemfl^"^ (see Introduction). The 
length ^ characterises the exponential decay of the eigenfunc¬ 
tions and is defined in terms of the re ciproc al of the Lyapunov 
exponent A. In ID A can be written aJSH]. 



'^R 



( 10 ) 


In Eq. the brackets denote averaging over both ensem¬ 
bles and the system unit cells, while the usual definition of the 
localization length considers only averages over ensembleJ^. 
The two definitions are equivalent. The relation between A 
and ^ is: 


ReA = 


( 11 ) 


where d is the mean length of the unit cell. The advantage of 
the approach based on the parameters Pn and Qn is that we 
can use polar (or action-angles) coordinates: 


Pn = Rn sin 0 


n? 


Qn = Rn COS 0 „ . 


( 12 ) 
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Without disorder, is a constant and 0^ increases by mi¬ 
nus the Bloch phase, —7, as we move from unit cell to unit 
cell. With disorder, the radius changes in every step, with 
Rn+i a function of Rn, ©n. and of the matrix elements of 
M^. The angle ©n+i only depends on 0^ and M^. For 
weak disorder a recurrence equation ( jB9| ) exists that, in the 
continuum limit, becomes a stochastic Ito equation which has 
a corresponding Fokker-Planck equatioiP^ . In this case, the 
first approximation for the density probability function of 0 ^ 
is uniform in the interval [ 0 , 27r] for 7 7 ^ 0 , 7 r/ 2 , tt. 

Writing Eq. ( p^ in terms of R and 0, and averaging over 
0 with uniform density probability, we obtain, up to second 
order md(j)i\ 

\=^(y^+Y 2 cos 2e„ + Fa sin 2e„ - - -^Y^^ , 

(13) 

where Yi, with i = 1,2,3 are defined in Appendix and 
depend on the matrix elements . 

In the following sections we will study the propagation of 
light through a disordered structure of alternating graphene 
sheets and dielectric layers. In this case each propagation ma¬ 
trix is determined by the widths Zi, the incidence angle Oi, 
the dielectric material parameters jii and Si, and the graphene 
conductivity a. In the present work we focus on the cases 
where disorder is present in the widths of the stacks (struc¬ 
tural disorder) and on graphene conductivities (compositional 
disorder) . Both are realistic situations that may occur in the 
fabrication of these structures. 

For the type of structural disorder studied here the width of 
each layer i of the cell is a random variable 


To proceed with the calculation of the Lyapunov exponent 
it is necessary to map the system parameters of the transfer 
matrix into the parametric matrix 0 .. There is not a 
unique way of doing this, but in what follows we make the 
simplest choice. 

1. Disordered photonic super lattice without graphene 

To model a disordered photonic super-lattice without the 
graphene layer we put / = 0 in Eq. and map (pi into the 
system parameters o^i, 0 ^ 2 , X, A (defined in Appendix [C|): 

sinh(/)i = sin 0^2, 

02 = 0(1 + arctan (x^ tan 0 ( 2 ), 
03=0(1+71/2, (16) 

where x =TE,TM. According to this mapping we can replace 
(pi in the expressions for Yi in Appendix and calculate the 
differentials dpi using Eq. © with Qi This will enable 

us to compute the Lyapunov exponent, given by Eq. <[13; the 
final result is 

A 2 

A = -^ (sin^ 0(2 kl af + sin^ ai cr|) , (17) 

2 sin 7 

which agrees with the result of Ref. 1^ for uncorrelated dis¬ 
order. The described procedure is repeated to calculate the 
Lyapunov exponents in the next sections. 

2. Disordered superlattice containing graphene layers 


Zi{n) = z^ + Ci{n), (14) 

where Q are uncorrelated random variables with zero mean 
and mean standard deviation ap ([Ci(^)]^) = o-f; z^ is the 
mean width of the i slab (the standard deviation should not 
be confused with graphene’s conductivity a). 

In the case of compositional disorder, the Eermi energy Ep 
is a random variable in each layer n 

EF{n)/h = LUpin) =ujpE Cp{n), (15) 

with (p Si random variable with zero mean and (C|^) = 

This determines how the graphene conductivity, given in Ap¬ 
pendix]^ is affected by disorder. 

In the next section we derive analytical expressions for A 
(Eq. in different regimes. To this end, we need to map pi 
in the system variables, calculate the differentials di, and use 
the results given in Appendix [B| 


A. Unit cell made of two different dielectric materials and a 
graphene sheet at the interface 

We consider a disordered superlattice composed of dielec¬ 
tric bilayers with a graphene sheet in between. The transfer 
matrix for the n unit cell is given by = {rri^i} and is 
explicitly derived in Appendix |C1[ 


The presence of graphene at the interface between the di¬ 
electrics results in a discontinuity in the tangential compo¬ 
nent of the magnetic field. The role of graphene on the op¬ 
tical properties of the superlattice increases as the value of the 
dimensionless parameter [3f f increases, with / = 
and given in Appendix We are interested in the loss¬ 
less regime in which the Bloch phase 7, given by Eq. ( |C3| ), is 
real. This regime sets up when (i) a (and therefore /) is a pure 
complex number and with i = 1 , 2 , is a pure real number; 
or (ii) cr is a pure real number so that evanescent propagation 
occurs in one of the layers. 

In the first case, we define = iB^ (see Appendix [^, 
where B is real, and we map the parameters pi in: 


sinh 01 = —B^ cos 0(2 + (A — D^) sin 0(2, 

02 = 0(1 + arg [A^ cos 0(2 + i(x + C'+) sin0(2] , 

03 = 0(1 + 7r/2 . (18) 


Eollowing the procedure of Sec. IIA1 the Lyapunov expo¬ 
nent is given by: 


A = 


2 sin^ 7 


{Kikfaf + Klklal) 


(19) 


where: 


Ki = -2/A*^|cosai + 


-A + 


sin 0(1, ( 20 ) 
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and / = if, K 2 is obtained by interchanging 1^2 and 
A ^ — A. Notice that if one plugs Eq. (|20|) with / = 0 into 
Eq. ( p^ , Eq. ^rf\ is obtained, as it should be. 


B. Unit cell made of one dielectric material and a graphene 
sheet at the interface 

Eor systems composed of bilayers of the same dielectric 
material with a graphene sheet in between, it is much easier 
to calculate the transfer matrix, which is given in Eq. ( |C6| ). In 
this case, the (j)i parameters read 

sinh(/)i = -A"r/, 

02 = 

03 = a + 7r/2, (21) 

Using Eq. ( p^ and the results of the Appendixj^we calculate 
the Lyapunov exponent for structures containing both random 
graphene conductivities (compositional disorder) and random 
widths (structural disorder), as detailed in the following. 


7. Compositional disorder 


Using the same procedure of subsection |II A 1[ we obtain 
the Lyapunov exponent: 



sin 2a 
sin 27 



2 


( 22 ) 


where ac is the fine structure constant and (jg is the mean 
standard deviation of the normalized graphene conductivity 


cr 


2 

9 


{a^) - (<7)2 
^2 


(23) 


2. Structural disorder 


Eor structural disorder where the stacks’ widths are given 
by Eq. the Lyapunov exponent reads 


A = 


2 sin^ 7 


(24) 


This concludes the analytical part of our work, which shall 
be compared to numerical simulations in the following sec¬ 
tion. 


III. NUMERICAL SIMULATIONS: RESULTS AND 
DISCUSSIONS 

A. Simulation procedure 

The numerical calculations are based on the transfer matrix 
method; the total transfer matrix for light propagating in a TV¬ 
layered system is 

(25) 


where the elements of are given by Eq. Transmis¬ 
sion is calculated by applying the boundary condition related 
to the fact that there is no incoming wave from the left: 


T = 



(26) 


and the localization length ^ is calculated by: 


f = -i(lnr), 


(27) 


where L = Nd and TV is the total number of unit cells with 
mean width d. The length L is chosen to be large enough 
to ensure the numerical calculation of the localization length 
converges. In the numerical procedure we first generate ran¬ 
dom variables Q (or (p) [see Eqs. ( p^ and ( p~5) )] from a uni¬ 
form distribution, and then calculate the transfer matrix using 
Eq. ( [^ . With the help of the results introduced in Appendix 
[C] we obtain the localization length using Eq. (^). The pro¬ 
cedure is repeated over nsampies and the mean value of the lo¬ 
calization length is calculated. We have verified that, for a 
sufficiently large TV, the value of ^ calculated for a single dis¬ 
order realisation coincides with its average over many disor¬ 
der realisations for smaller systems; in other words, we have 
verified that ^ is a self-averaging quantity. Eurther details of 
the transfer matrix method are given in Appendix 


B. Results 


Light transmission depends on the graphene conductivity a 
and on the medium impedances, defined as ( see Appendixj^: 

•te _ ^TM ^ VMi cos Oi. (28) 




l^i 


We shall focus in the lossless regime with cos 7 = 0 
and 5Recos7 < 1. Erom Eq. (C3), this regime occurs when¬ 
ever / (and consequently a) is a pure complex number or for 
= 0, in which case one of the slabs supports a evanescent 
mode. When the Drude term dominates, the imaginary part of 
the conductivity is positive (see Appendix]^. Eor frequencies 
slightly below 2 ujf, the inter-band term dominates and the 
imaginary part of the conductivity is negative (see Appendix 
0- When the frequency becomes larger than 2ujf, the imagi¬ 
nary part goes to zero and the real part tends to ao = e^/4/i. 

In the following numerical calculations, random variables 
have a uniform distribution with (x C [—Tx/2^Tx/2], with 
X = 1, 2 for structural disorder and x = F for compositional 
disorder. 


C. Drude regime when: ^ea ^ 0, ^ma > 0 

When ujpf^ ^ ^ ^F (where T is the broadening en¬ 

tering in the conductivity), graphene conductivity can be ap¬ 
proximated by (see Appendix |D1| ): 

4 ujf 

a = Ido -. 

TT UJ 


M = n7iM". 


(29) 
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For Ef ~ 0.3 eV (a typical value for the graphene Fermi 
energy), the range of frequencies corresponds to the infrared 
spectral regime. In the following we focus in three regimes: 
impedance matching in the double layered system [Zi = Z 2 , 
in Eq. ( [28] )] with structural disorder, compositional disorder 
in one layered system, and the attenuated field regime (ATR) 
with structural disorder. 


1. Impedance matching in two-layered system with structural 
disorder 


Using the Snell-Descartes law, Eq. ( |C2| ), and the 
impedances in Appendix one can verify that for materi¬ 
als without magnetic response (/ii = /i 2 = 1), there is no TE 
mode that allows the impedance matching. In the TM mode 
the impedance matching occurs when the angle of incidence in 
layer 1 obeys the relation sin^ 0 i = ^ 2/(^1 +^ 2 ), ei ^ 62 . 

When Zi = Z, /3i = /3, it follows from Eqs. ( [T^ and Eq. 
( | 20 | ) that: 


A = 2 


TTC sin 7 


Sifii cos^ ai cos^ 

i=l 


(30) 


where we neglected the term P in comparison to / (which 
in the Drude regime is always valid for a sufficient large uj). 
In this case, in Eig. [^ the localization length ^ is calculated, 
both analytically and numerically, as a function of frequency. 
The dispersion relation is also shown in Eig. |^. It is impor¬ 
tant to point out that the agreement between the analytical and 
numerical calculations is very good, except when 7 approach 
0 or TT. This is due to the fact that, in the analytical deriva¬ 
tion of the Lyapunov exponent, the recurrence equation ( |B9| ) 
is ill defined at these points, so that the distribution of random 
variables is not uniform. Remarkably, Eig. [^ reveals that in 
the impedance matching regime, ^ does not follow the well- 
known asymptotic power law behaviour for low frequen¬ 
cies. Rather, ^ exhibits a periodic dependence on uj for low 
frequencies, a result that is intrinsically related to the graphene 
conductivity properties. Indeed, it can be explained by the fact 
that the linear increase of the wavenumber with frequency is 
cancelled by the simultaneous decrease of graphene’s conduc¬ 
tivity (Drude term, see Eq. [^, which scales with 1/uj. The 
periodicity in ^ follows from the periodicity in the dispersion 
relation, shown in Eig. |^. Eor the lossy and Drude regimes, ^ 
approaches the same value as the frequency increases, and the 
real part of the Drude conductivity goes to zero. 

Eigurej^ shows ^ as a function of frequency for two differ¬ 
ent values of the incidence angle 0. It reveals that the presence 
of graphene layers has also an important effect in the so-called 
Brewster modes in disordered systems. In ID disordered op¬ 
tical systems, the so-called Brewster modes occur at some 
specific frequencies and incident angles for which j rey hes 
anomalously high values, larger than the system size^^EH, por 
non-magnetic (/ii = /i 2 = 1 ) superlattices made of positive 
refractive-index media, these anomalously delocalised modes 
arise from the suppression of reflexion at the interface s of a ID 
disordered system illuminated by a TM incident wavJ^^EH^ As 


a result, the system becomes fully transparent. The presence 
of graphene induces additional reflections at each interface of 
the superlattice, resulting in an attenuation of this Brewster 
mode, as it can be seen from fig. [^ 




FIG. 2. Color on-line, (a) Dispersion relation in the impedance 
matching regime and TM mode, (b) Localization length as a func¬ 
tion of frequency with Ti = 5/im, = 1.2mm, Ep = 0.2 eV, 

r = 260./ieV, N = 5000, risampies = 100, si — gn — gi 2 — 1, 
82 — ‘^,0 — 7r/3.The solid line in (b) is the analytical result, whereas 
the dots correspond to two different numerical simulations for differ¬ 
ent regimes of the optical conductivity of graphene: (i) a — ^maD 
(red points) and (ii) a = an E cfi (blue points). 


2. ATR regime in one-layered system 

The plasmon-polariton mode in graphene can be excited for 
example, by a prism in the Otto configuration 1^ . This is the 
regime we will explore in this section. We consider a periodic 
array of graphene/air unit cells (medium 2 ) in between a di¬ 
electric (medium 1). In this case the total transfer matrix M 
is obtained considering the boundaries between the prism and 
the superlattice: 

M = Mi^2Yl{Mj)M2^i, (31) 

3 

where Mi ^2 refers to the transfer matrix describing light 
propagation from the medium 1 (dielectric) to medium 2 (air); 
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FIG. 3. Color on-line. Localization length as a function of frequency 
in the impedance matched regime for two values of incidence angle: 
0 = 60° (blue line), 0 = 45° (red circles). = 0.5/im, = 

120/im, Ef — 0.2 eV, V — O./xeV, N — 5000, risampies = 100. 



FIG. 4. Color on-line. Localization length as a function of incidence 
angle in the impedance matched regime at the vicinities of a Brewster 
mode. Ti — 50/im, — 120/im, N — 5000, risampies = 100. 

The black solid line corresponds to the grapheneless case; red circles 
correspond to the case where graphene is present in the superlattice 
{Ef = 0.2eVandr = OeV.) 


M 2 ^i refers to the reverse propagation. Mj is the trans¬ 
fer matrix of the unit cell air/graphene with random widths 
(medium 2). 

From Eq. one can see that for the evanescent mode a 
is a pure complex number and the first term in the right hand 
side becomes a hyperbolic cosine, which is greater than 1 for 


any a. As a. result, the Bloch phase is real only if the second 
term in the right hand side of |C7| is negative. This situation 
occurs for pure positive complex /; in this case (3 is also a 
pure positive complex number, which is only possible in the 
TM mode [see Eqs. and (|C4|)]. 

Eor an incidence angle 0i above the critical angle for total 
reflection at the interface 1/2, a plasmon-polariton can be ex¬ 
cited, allowing for frustrated total internal reflection. In this 
case light propagation occurs due to the presence of periodic 
graphene sheets. The effective impedance in the medium 2 
depends on the properties of the layer 1 as: 



In Eig. [^the localization length is calculated in the ATR 
regime using both numerical and analytical methods; the 
agreement is excellent. In the Drude regime ^ is inversely pro¬ 
portional to the Eermi energy. Also shown is the localization 
length when the dielectric necessary to excite the ATR field 
is removed; we call this situation the normal field. The ATR 
field is characterized by exponentials with argument Eujnzjc. 
When the frequency increases and the length c/ tiuj becomes 
smaller than the width z of the dielectric slab (air in this case) 
light propagation comes to a halt, as the plasmon-polariton 
localized in a graphene layer cannot excite the adjacent layer. 
We can see that the ATR for the parameters of Eig. fills the 
band gap of the normal field. Also the increase of disorder 
implies in the decrease of as expected. 

Notice that ignoring the interband term and making Ep = 0 
is equivalent to remove the graphene sheets, therefore making 
disorder in random widths of air meaningless. Hence the lo¬ 
calization length diverges, as can be seen in Eq. ( [24| ), where 
cr 0 implies in a vanishing Lyapunov exponent. 

3. One layer system with compositional disorder 

In the compositional disorder regime and for the one lay¬ 
ered system, ^ decreases as /3 increases. Eor the TE mode, 
P can only be greater than 1 for materials with magnetic re¬ 
sponse, /i > 1. Eor the TM mode, p is proportional to the 
dielectric constant and to cosO, thus for grazing incidence, 
the system becomes fully opaque. 

In the Drude regime the asymptotic behaviour of the lo¬ 
calization length goes as This can be understood as fol¬ 
lows: as the frequency increases the graphene conductivity 
decreases as and thus the influence of the graphene layer 
disappears. 

The effect of compositional disorder is shown in Eig. 
where the Eermi energy is randomly distributed around the 
mean value E^ = 0.6 eV. ^ is inversely proportional to the 
mean standard deviation of the Eermi energy. We study the 
effect of increasing absorption in graphene layers, which de¬ 
pends on the real part of the conductivity and is proportional 
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FIG. 5. Color on-line, (a) Dispersion relation for the ATR regime 
(red) and for normal field (blue), (b) Localization length as a function 
of frequency for = 12/im, Ep = 0.1 eV, 0i = tt/S, si = 
2,/ii = = 1, iV = 50000, risampies = 1- Thc yellow 

circles (green squares) and orange diamonds (blue triangles) refer to 
the ATR (normal) field with T = 0.5/imandT = 5/im, respectively. 
The cyan and purple lines refer to the analytical approximation. 


to the relaxation rate F. The length ^ decays rapidly when 
the frequency reaches 2 ujf, and interband transitions start to 
occur, an effect that may be related either to absorption or 
to Anderson localization. The numerical calculation is per¬ 
formed with the full graphene conductivity (Drude plus inter¬ 
band) and then compared to the case where only the Drude 
term is present. The analytical approximation is calculated 
with the Drude term only, and agrees very well with the nu¬ 
merical simulation except at the band edges 7 = 0,7r/2,7r. As 
already discussed, this disagreement is related to the fact that 
the probability distribution of 0^ is not uniform for these val¬ 
ues of 7 . The analytical approximation has a peak at 7 = 7r/2 
(see denominator of Eq. [2^ . The numerical calculations show 
that near the band gap (7 = 0, tt see Eq. ^ goes to zero, 
and the peak at 7 = 7r/2 does not occur. 




FIG. 6. Color on-line, (a) Real and imaginary parts of the graphene 
optical conductivity in the compositional disordered case, a = aoE 
CT/, and the Drude conductivity an when F = 0. (b) Localization 
length as a function of frequency with = 1.2/im, 0 = tt/A, e = 
= 1, Ep = 0.6 eV, Tf = 0.12 eV, N = 5000 and increasing 
relaxation rate F. up ~ 3 x 10^^ Hz. The blue triangles (and solid 
blue line) refer to a calculation where only the Drude conductivity 
with F = 0 is used. The other data sets refer to the use of the full 
optical conductivity of graphene with different F values. 


D. Complex interband regime when: ^ea ^ 0, ^ma < 0 

When uo < 2 ccf, the imaginary part of the optical conduc¬ 
tivity of graphene becomes negative and can be approximated 
by 


a = Id T 


1(JQ- 


4 (jjp 


(34) 


where aj is given by Eq. ( |D4[ ). In this case the imaginary part 


of (j becomes negative, and the ratio between the imaginary 
and real parts of cr becomes lower than in the Drude regime 
for typical values of F and Ep. Therefore, in this case the ex¬ 
ponential decay of transmission is essentially due to absorp- 
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tion rather than to Anderson localization. Therefore, in this 
case, our approach for studying Anderson localization using 
the localization length is inadequate. It is worth commenting 
that experimentally it is possible to distinguish between ab¬ 
sorption and Anderson localization by investigating the vari¬ 
ance of the normalized total transmission, as proposed in Ref. 
Ell. For a one-layered system in the ATR regime with trans¬ 
fer matrix given by Eq. (C6), the change in the sign of / has 
qualitatively the same effect in the dispersion relation ( |C7| ) of 
interchanging TE and TM modes, which changes the sign of 
/3. 


When the frequency becomes larger than 2uof, the real 
part of the conductivity approachs do while the imaginary 
part vanishs. In this regime, the role of the graphene sheets 
consists, essentially, in absorbing light leading to a vanishing 
transmission after few stacks. 


IV. CONCLUSIONS 

In conclusion, we have investigated light propagation in 
ID disordered superlattices composed of dielectric stacks and 
graphene sheets in between. We introduced disorder either 
in the graphene material parameters (compositional disorder), 
such as the Fermi energy, or in the widths of the dielectric 
stacks (structural disorder). For both cases we derived an an¬ 
alytical expression for the localization length ^ and compared 
the results with numerical calculations based on the transfer 
matrix method. A very good agreement between numerics and 
the analytical expression was found. We demonstrated that, 
for structural disorder and when the impedances of the layers 
are equal, the localization length does not follow the well- 
known asymptotic behaviour ^ (X UJ Rather, it exhibits an 
oscillatory dependence on frequency, as a result of the pres¬ 
ence of the Drude term in the graphene conductivity. Also in 
the impedance matching regime, we show that graphene has 
an important impact on the Brewster modes, anomalously de¬ 
localised modes at given frequencies and incident angles at 
which ^ diverges. Indeed, the presence of graphene induces 
additional reflections inside the disordered medium, leading 
to a strong attenuation of the Brewster modes. We investi¬ 
gated how intra and interband transitions in the graphene con¬ 
ductivity impact on identifying the regimes where Ander¬ 
son localization and absorption dominates light transmission. 
Altogether, our findings unveil the role of graphene on An¬ 
derson localization of light, paving the way for the design of 
graphene-based, disordered photonic devices in the THz spec¬ 
tral range. 
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Appendix A: Matrix transformation 

The relation can be interpreted as a dis¬ 

crete set of points in the phase space 'ipL- With the trans¬ 
formation Mreafi 


-^real 


1 

2 


1 — i 
-1 + 


i 



(Al) 


the matrix is now real, and defining = 

we have in the phase space that in the system 

without disorder the trajectory is given by a ellipse. Prom this 
we can find a transformation Mcirde to a circle: 


Afcircle — 


COST 


V smr 


1-1 


—V "smr rcosT 


(A2) 


where: 


V 


2 


sin 7 

cosh (j)\ sin 02 + sinh (j)\ ’ 



and making [Q P]^ = MdrcieV^', 


(A3) 


fQn\ _ ( ^cosT rsinr\ fxn 

\Pn) ~ l^-r“^sinT rcosry \yn 


(A4) 


Appendix B: Lyapunov Exponent 


The Lyapunov exponent is given by: 



where 



(Bl) 


1^1 = 


1 

sin^ 7 


with: 


(B2) 


f/i = 2sin^02, (B3) 

1/2 = 2 sinh^ 05 cos^ 7, (B4) 

1/3 = 2 sinh^ 05 sin^ 7, (B5) 

f/4 = — sinh205 sin205, (B6) 


I2 = [—2 sin 0^(501 + cos 0^ sinh 205 ((502 — (503)] , (B7) 
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sinh (pi (cos^ 7(502 + sin^ 7 ( 503 ) — cos 7 sin 02 ( 50 i 

^3 = 2-^- 

— sin 7 

the angle 0 obeys the recurrence equation: 

0n+l = 0n - 7 + CSC 7, 

with: 

Cn = [cos 7 sinh 0^(502 — sin 02 ^ 0 i] cos ( 20 ^ ~ 7 ) + 

sinh 05 cos 7 sin(20^ — j)S(l)s , (BIO) 


(B8) 

(B9) 


Appendix C: Photonic Crystal 

1. Unit cell made of two different dielectrics and a graphene 
sheet at the interfaces 


with A™ = +1,AT^ = -1. 


2. Unit cell made of one dielectric and a graphene sheet at the 
interface 


When there is only one dielectric, with width z and 5 , fi per- 
missivity and permeability, intercalated by graphene sheets, 
the transfer matrix is given by: 


M = 




(1 + p^f)e-^^) ’ 


(C6) 


The transfer matrix whose elements arfl^ 

m{i = [A^ cos a2 + i(x + sin 0^2] , 

^12 = cos a 2 i {A D^) sin a 2 ] , 

= [—B^ cos a 2 — i {A D^) sin a 2 ] , 

m ^22 ~ [^+cos (3^2 —(x + C'-) sina 2 ] (Cl) 

where x = TE, TM and the diverse parameters are given in 
appendix 

The Snell-Decartes law hold: 

y/eifii sin 0i = ^/S 2 |J ^2 sin O 2 , (C2) 

and the dispersion relation is given by: 

cos 7 = cos (3^1 COS (3^2 ~ (x + sin Q^i sina2 

+2i/ (01 cos ai sin 0^2 + 0f cos 0^2 sin ai)(C 3 ) 


where a = y/JIez cos 0 with the dispersion relation: 

cos 7 = cos a — ij3^ f sin a. (C7) 


Appendix D: Graphene Optical Conductivity 

For completeness we give here the expressions for the op¬ 
tical con ductiv ity of graphene, whose derivation can be found 
elsewherd^^BU xhe graphene optical conductivity of graphene 
is a sum of a Drude term, cfd, and an inter-band contribution, 
a I, reading: 

(J = cfd^cfi, (D1) 


where: 


=zl^, = 


ZTE^ 


ki = .ysiiiiUi/c cos di, 

OCi 

Al = {l±2fp^), 

= 2/A“/?f, 

Cl = ±2//3f + 

= 2fX^p^p^, 


A' = 5 (•?“ - 


X = 


/ = 


ac/iQ 


zr = ^^cosei 


z7^ = 


l^i 


COS dj 


(C4) 


(C5) 


where the Drude term is given by: 

(7 D 4cJ F 1 


(Jo TT r — icC ’ 

and the interband term aj = a'j icf'l have the real part 


(D2) 


= 1 


O-Q 


1 uj — 2 ujf 1 ^ 2ujf 

- arctan---arctan--- 

r TT r 


TT 


and the imaginary part 


(jj _ ^ In ^ ^ 


2tt {2ujf — + T^ ’ 


(D3) 


(D4) 


with the Fermi frequency given by 


ujf — 


\Ef\ 


(D5) 
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